import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.preprocessing import scale
from sklearn.impute import SimpleImputer
pd.set_option("display.max_columns",500)
pd.options.mode.chained_assignment = None  # default='warn'

def newHpFilter(path):
    initExcel = pd.read_excel(path, sheet_name='initial', header=None)
    # 数值部分：[5:,1:]
    initData = initExcel.loc[5:, 1:]  # initData的行标是和initExcel一样的，但是列标重新从1到76编号
    hpT = pd.DataFrame()  # 用来装载hp滤波后的dataframe
    subMat=pd.DataFrame() #装cycle的
    notNa = ~initData.isna()
    for col in range(1, initData.shape[1] + 1):
        cycle, trend = sm.tsa.filters.hpfilter(initData.loc[notNa.loc[:, col], col],lamb=1600)
        subMat=pd.concat([subMat,cycle],axis=1,ignore_index=True)
        hpT = pd.concat([hpT, trend], axis=1,ignore_index=True)
    #return subMat #行标：5:76 列标：1：2159
    stdMat=subMat.std(ddof=0,axis=0)
    newSubMat=subMat.copy(deep=True)
    for ncol in range(1,subMat.shape[1]+1):
        for nrow in range(5,subMat.shape[0]+5):
            if subMat.loc[nrow,ncol]>1.5*np.abs(stdMat.loc[ncol]):
                newSubMat.loc[nrow, ncol] = 1.5 * np.abs(stdMat.loc[ncol])
            elif subMat.loc[nrow,ncol]<-1.5*np.abs(stdMat.loc[1]):
                newSubMat.loc[nrow, ncol] = -1.5 * np.abs(stdMat.loc[ncol])

    # newSubMat.to_excel(excel_writer='D:\\subMat.xlsx')
    # hpT.to_excel(excel_writer='D:\\hpT.xlsx')
    resMat=newSubMat+hpT #resMat 任然是index:[5,76]  columns:[1,2159]
    fAns=pd.DataFrame()
    for rCol in range(1, initData.shape[1] + 1):
        Cycle1, Trend1 = sm.tsa.filters.hpfilter(resMat.loc[notNa.loc[:, rCol], rCol], 10)
        fAns=pd.concat([fAns,Trend1],axis=1)
    #fAns.to_excel(excel_writer='D:\\shuju.xlsx',sheet_name='sheet1',header=None,index=False)
    initExcel.loc[5:, 1:] = fAns
    # 交换第0行和第4行位置，让X1,X2,X3...这一行变成第一行
    a, b = initExcel.loc[0, :].copy(deep=True), initExcel.loc[4, :].copy(deep=True)
    initExcel.loc[0, :], initExcel.loc[4, :] = b, a
    # 处理日期格式，不然写入excel后日期就变成了#######
    initExcel.loc[5:, 0] = initExcel.loc[5:, 0].apply(lambda x: str(pd.to_datetime(x).date()))
    # 指定header=None，Index=False，不然会将行、列的索引号也写入excel中
    initExcel.to_excel(excel_writer='D:\\SoftWareHomeWork2\\in新的python跑出来hp滤波.xlsx', sheet_name='Sheet1', header=None, index=False)
    return initExcel


#cross_cor
# 测试代码：
# x=pd.DataFrame(np.array([3.98,8.23,9.9,3.3,2.95,21.43,17.64,8.222]))
# y=pd.DataFrame(np.array([1.38,2.90,4.34,3.64,5.89,7.56,8.23,9.53]))
# L=6
# ansT=cross_cor(x,y,L)
# print(ansT)
#经过和matlab的多组不同数据结果对比，该函数正确
def cross_cor(x,y,L)->pd.DataFrame:#x,y都是Series
    T=y.shape[0]
    p=pd.DataFrame(np.zeros((1,2*L+1),dtype=np.float))#1行，2L+1列
    varx=(x-x.mean())*(x-x.mean())
    # sigma_x=varx.sum().apply(lambda x:x**0.5)[0]#sum()返回的是一个1*1的df，apply后还是1*1
    sigma_x= np.sqrt(varx.sum())#sum()返回的是一个1*1的df，apply后还是1*1
    vary=(y-y.mean())*(y-y.mean())
    # sigma_y = vary.sum().apply(lambda y: y ** 0.5)[0]
    sigma_y = np.sqrt(vary.sum())
    print(sigma_x)
    print(sigma_y)
    j=0
    for i in range(-L,0):#取不到0
        n1=T+i
        n2=-i
        u=(y.iloc[0:n1,:]-y.mean())
        u=u.values
        v=x.iloc[n2:,:]-x.mean()
        v=v.values
        uv=u*v
        # m1=((y.iloc[0:n1,:]-y.mean())*(x.iloc[n2:,:]-x.mean()))#这样写不行，因为(y.iloc[0:n1,:]-y.mean())和x.iloc[n2:,:]-x.mean()的index不同，*会变成Nan
        s1=uv.sum(axis=0)
        # ss1=s1[0]
        p[j]=pd.DataFrame([s1[0]/(sigma_y*sigma_x)])
        # print(p[j])
        j+=1
    for i in range(0,L+1):
        n1=T-i
        n2=i
        e=(x.iloc[0:n1,:]-x.mean())
        e = e.values
        r=(y.iloc[n2:,:]-y.mean())
        r=r.values
        er=e*r
        s2=er.sum(axis=0)
        # print(s2[0])
        p[j]=pd.DataFrame([s2[0]/(sigma_x*sigma_y)])
        # print(p[j])
        j+=1
    period=pd.DataFrame(np.arange(L,-L-1,-1)).T
    p=pd.concat([period,p],axis=0,ignore_index=True)#纵向拼接
    p=p.T
    p=pd.concat([p,p[1].abs()],axis=1,ignore_index=True)#增加一列abs值，方便查找至最大的那一列，待会不取就完事
    #现在的p第0列是期数，第1列是对应的相关系数，第2列是第1列的绝对值，要找第2列的最大值的索引
    # print(p)
    #新的
    # p=p.fillna(value=np.NINF)
    p = p.values
    idx = p[:, 2].argmax(axis=0)
    # print(idx)
    ansT = p[idx, [1, 0]]
    ansDf = pd.DataFrame(ansT)
    ansDf = ansDf.T
    # print(ansDf)
    return ansDf

#求解相关系数
def getCorre(path,L):
    excel=pd.read_excel(path,sheet_name='Sheet1',header=None)#header=None:代表不使用源表的表头
    numData=excel.iloc[5:,1:]
    numData=numData.reset_index(drop=True)#注意，现在的numData行索引从0开始，列索引从1开始，因为没有重置列索引
    print(type(numData))
    [h,k]=numData.shape
    rho=pd.DataFrame(np.zeros((k,2))*np.nan)
    xnan=~numData.isnull()
    for col in range(1,k+1):#[1,k]，能取到k的
        #nx的index还是原来的，即如果非空的行为12,13...，那么nx的第一行的index就是12
        nx=numData.loc[xnan[col],col].values#iloc是按照索引来的，索引[:,1]取的是第2列而不是第一列，而loc[:,1]就是取得第一列
        ny=numData.loc[xnan[col],1].values#注意这里ny也要取xnan[col]，这样才是对应的
        nx=(nx-np.mean(nx))/np.std(nx,axis=0)
        ny=(ny-np.mean(nx))/np.std(ny,axis=0)
        x=pd.DataFrame(nx)
        y=pd.DataFrame(ny)
        res=cross_cor(x,y,L)
        # print(res)
        res=res.values
        res=pd.DataFrame(res).loc[0]
        rho.loc[col,:]=res
    rho_sorted=rho.sort_values(by=[1,0],axis=0)
    flag=rho_sorted.index
    name=excel.iloc[0:5,1:].T
    name_sorted=name.loc[flag,:]
    ansDf=pd.concat([name_sorted,rho_sorted],axis=1,ignore_index=True)
    ansDf.to_excel(excel_writer='D:\\SoftWareHomeWork2\\new求解相关系数.xlsx',sheet_name='Sheet1',header=None,index=False)
    return ansDf
    #测试数据
    # path = 'D:\\pyhpfilter.xlsx'
    # L = 6
    # ansDf = getCorre(path, L)
    # print(ansDf)

#获取每一期绝对值大于0.6且最大的那num个
def getMaxAbs(path,num,L):#path是文件路劲，num是取前num个，L是期
    # path='D:\\ProgramData\\MATLABProjects\\'
    table=pd.read_excel(path,sheet_name='Sheet1',header=None)
    table=table.rename(columns={0:'X标号',1:'指标名称',2:'频率',3:'单位',4:'指标ID',5:'相关系数',6:'期'})
    table['sort']=table['相关系数'].abs()
    table=table.sort_values(by=['期','sort'],ascending=[True,False])
    table=table.drop(labels=['sort'],axis=1)
    table.to_excel(excel_writer='myTestBest.xlsx',header=True,sheet_name="Sheet1")

    ans=pd.DataFrame(columns=['X标号','指标名称','频率','单位','指标ID','相关系数','期'])
    for i in range(-L,L+1):
        tmp=table.loc[(table['期']==i) & (table['相关系数'].abs()>=0.6)]
        # tmp=table[table['相关系数']>=0.6 and (table['期']==i)]
        if tmp.shape[0]>num:
            tmpAns=tmp.head(num)
            ans=pd.concat([ans,tmpAns],axis=0,ignore_index=True)
        else:
            ans=pd.concat([ans,tmp],axis=0,ignore_index=True)
    print(ans)
    ans.to_excel(excel_writer='D:\\SoftWareHomeWork2\\选取相关度最高的888个.xlsx',header=None,sheet_name='Sheet1',index=None)

def myPCA(path1,path2):#path1是hp滤波后的，path2是相关系数表
    excel1=pd.read_excel(io=path1,sheet_name='Sheet1')#以第一行作为表头
    excel2=pd.read_excel(io=path2,sheet_name='Sheet1',header=None)
    excel2=excel2.rename(columns={0:'X标号',1:'频率',2:'单位',3:'指标ID',4:'指标名称',5:'相关系数',6:'期'})
    group=excel2.groupby('期')#分组聚集，size()方法返回每组的行数
    gSize=group.size()#gSize是一个Series
    gNum=gSize.shape[0]#gNum是组数，等于2*L+1
    L=(gNum-1)//2
    pos=0
    writer = pd.ExcelWriter('D:\\SoftWareHomeWork2\\fiannly.xlsx')
    for t in range(0,gNum):
        #每组数据的行标：[pos,pos+gSize.iloc[t]),每次执行完后pos=pos+gSize.iloc[t]
        Xids=excel2.iloc[pos:pos+gSize.iloc[t],0]
        XidsT=Xids.T#转置一波，方便直接在excel1里面利用列名查找
        goalCols=excel1[XidsT]#这就是我们要找的这一期的那些列
        for row in range(pos,pos+gSize.iloc[t]):
            if excel2.iloc[row,5]<0:#这个变量的相关系数小于0，那么那一列数据去相反数
                Xid=excel2.iloc[row,0]
                goalCols.loc[4:,Xid]=-goalCols.loc[4:,Xid]
        goalData=goalCols.iloc[4:,:]#取数值部分,index和goalCols一样，从4开始，columns变成指定的列
        notNaMat=goalData.isna()
        goalDataNp=goalData.values#转出numpy数组
        imp_mean = SimpleImputer(missing_values=np.nan, strategy='mean')
        imp_mean.fit(goalDataNp)
        goalDataNp=imp_mean.transform(goalDataNp)
        npGoalData=scale(goalDataNp)#标准化，标准化之后就自动根据协方差矩阵进行主成分分析了
        pca = PCA()  # 可以调整主成分个数，n_components = 1
        pca.fit(npGoalData)
        # print(pca.explained_variance_)  # 输出特征根
        # print(pca.explained_variance_ratio_)  # 输出解释方差比
        # print(pca.components_)  # 输出主成分
        pcaComponent=pca.components_#主成分矩阵
        pcaFirstCol=pcaComponent[:,0]#，只取第一列
        pcaFirstRow=pcaFirstCol.T
        ansNp=pcaFirstRow*goalDataNp
        ansNp=ansNp
        resDf=pd.DataFrame(ansNp)#此时resDf的index已经从0开始了，原来的是从4开始，双射，逆映射

        #开始处理文本部分为正确格式
        text=goalCols.iloc[0:4,:]#文本部分
        text=text.T#通过转置删除原来的columns，不然拼接不成功
        text=text.reset_index(drop=True)
        text=text.T
        header=goalCols.columns
        header=pd.DataFrame(header).T
        ansText=pd.concat([header,text],ignore_index=True)
        #至此，成功得到了正确格式的文本部分
        #拼接文本和数值部分，两者的列索引是一样的，所以不用改,直接拼接，然后写入文件
        ans=pd.concat([ansText,resDf],axis=0,ignore_index=True)
        gdpColName=excel1.columns.tolist()[1]
        other=excel1.iloc[0:,1]
        gdpColName=pd.DataFrame([gdpColName])
        gdpCol=pd.concat([gdpColName,other],axis=0,ignore_index=True)
        ans=pd.concat([gdpCol,ans],axis=1,ignore_index=True)
        ans.to_excel(writer,sheet_name='Sheet'+str(t-L),header=True,index=None)
        print('pos='+str(pos))
        pos = pos + gSize.iloc[t]
    writer.save()
    writer.close()
    return

if __name__ == '__main__':
    path='D:\\2020data01.xlsx'
    newHpFilter(path)
    getCorre('D:\\SoftWareHomeWork2\\in新的python跑出来hp滤波.xlsx',4)
    getMaxAbs('D:\\SoftWareHomeWork2\\new求解相关系数.xlsx',num=10,L=4)
    myPCA('D:\\SoftWareHomeWork2\\in新的python跑出来hp滤波.xlsx','D:\\SoftWareHomeWork2\\选取相关度最高的888个.xlsx')
    print('跑完了')
    # path = 'D:\\in新的python跑出来hp滤波.xlsx'
    # L = 6
    # ansDf = getCorre(path, L)
    # x=pd.DataFrame(np.array([np.nan ,5.3, 6.9 ,8.23 ,4.24,6.64,8.43,9.32,5.65,3.112]))
    # y=pd.DataFrame(np.array([28,  25,  7 , 0,  24,  6 , 23 , 16,  3 , 3]))
    # L=6
    # ansDf=cross_cor(x,y,L)
    # print(ansDf)
    # p1='D:\\rightHp.xlsx'
    # p2='选取相关度最高的888个.xlsx'
    # myPCA(p1,p2)
    # path='D:\\new求解相关系数.xlsx'
    # getMaxAbs(path,10,6)